# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to run main models with religion control and create:
#   - Table A.6 (logit results with religion control)

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data
myData <- readRDS("data/afrobarometer.rds")

# Explanatory variables:
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)

# MODELS

# Model 1: Media Aggregate (Controls + CFE + DCSE)
model <- 
  sexuality2 ~
  media_aggregate +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.1.all <- glm(model, family = binomial, data = mdata)
logit.result.1.all.robust <- coeftest(logit.result.1.all, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.1.all.robust
lm.result.1.all <- lm(model, data = mdata)
lm.result.1.all.robust <- coeftest(lm.result.1.all, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.1.all.robust

# Model 2: Media Types (Control other media + Controls + CFE + DCSE)

# Radio
model <-
  sexuality2 ~
  radio + 
  media_noradio +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.2.radio <- glm(model, family = binomial, data = mdata)
logit.result.2.radio.robust <- coeftest(logit.result.2.radio, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.2.radio.robust
lm.result.2.radio <- lm(model, data = mdata)
lm.result.2.radio.robust <- coeftest(lm.result.2.radio, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.radio.robust

# TV
model <-
  sexuality2 ~
  tv + 
  media_notv +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.2.tv <- glm(model, family = binomial, data = mdata)
logit.result.2.tv.robust <- coeftest(logit.result.2.tv, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.2.tv.robust
lm.result.2.tv <- lm(model, data = mdata)
lm.result.2.tv.robust <- coeftest(lm.result.2.tv, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.tv.robust

# Newspaper
model <-
  sexuality2 ~
  newspaper +
  media_nopaper + 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.2.newspaper <- glm(model, family = binomial, data = mdata)
logit.result.2.newspaper.robust <- coeftest(logit.result.2.newspaper, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.2.newspaper.robust
lm.result.2.newspaper <- lm(model, data = mdata)
lm.result.2.newspaper.robust <- coeftest(lm.result.2.newspaper, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.newspaper.robust

# Internet
model <-
  sexuality2 ~
  internet +
  media_nointernet +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.2.internet <- glm(model, family = binomial, data = mdata)
summary(logit.result.2.internet)
logit.result.2.internet.robust <- coeftest(logit.result.2.internet, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.internet <- lm(model, data = mdata)
lm.result.2.internet.robust <- coeftest(lm.result.2.internet, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.internet.robust

# Social Media
model <-
  sexuality2 ~
  social_media +
  media_nosocialmedia +
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  as.factor(country) + 
  as.factor(religion)
mdata <- extractdata(model, myData, na.rm=TRUE, #only extract data used in model
                     extra = ~respno + district) #add respno & district
logit.result.2.social_media <- glm(model, family = binomial, data = mdata)
logit.result.2.social_media.robust <- coeftest(logit.result.2.social_media, vcov. = function(x) cluster.vcov(x, ~ district))
logit.result.2.social_media.robust
lm.result.2.social_media <- lm(model, data = mdata)
lm.result.2.social_media.robust <- coeftest(lm.result.2.social_media, vcov. = function(x) cluster.vcov(x, ~ district))
lm.result.2.social_media.robust


# Generate Latex table with logit results: Table A.6 in Appendix 
stargazer(logit.result.1.all.robust, logit.result.2.radio.robust, logit.result.2.tv.robust,
          logit.result.2.newspaper.robust, logit.result.2.internet.robust, logit.result.2.social_media.robust,
          type = "latex",
          no.space=TRUE, 
          dep.var.caption = "DV: Homosexual as Neighbor (0: dislike, 1: don't care/like)",
          dep.var.labels.include = FALSE,
          covariate.labels = c("Media aggregate", "Radio", NA, "TV", NA, "Newspaper",
                               NA, "Internet", NA, "Social media", NA, "Female", 
                               "Education", "Religiosity", "Age", "Income", "Urban"),
          title = "Effect of Media Consumption on LGBT Attitudes (Logit Models)",
          omit = c("country", "religion"),
          keep.stat = c("n", "AIC"),
          notes = c("All models include country fixed effects. Standard errors are clustered at the district level.")
)
